package model.util;

import java.awt.image.BufferedImage;
import java.util.ArrayList;

import model.doe.AbstractDoe;

import org.apache.commons.math3.transform.DftNormalization;
import org.apache.commons.math3.transform.FastFourierTransformer;
import org.apache.commons.math3.transform.TransformType;

import test.Complex;

import edu.emory.mathcs.jtransforms.fft.DoubleFFT_2D;

public class FourerExecutor {

	public enum DoeType {
		AMPLITUDE, PHASE, PHASE_MODAN_AMPLITUDE
	}
	
	
	double[][] amplMatrix = null;
	double[][] phaseMatrix = null;
	double[][] resultAmpl = null;
	double[][] resultPhase = null;
	
	
	/**
	 * 
	 * @param amplMatrix - amplitude matrix for fourer transform (module)
	 * @param phaseMatrix - phase matrix for fourer transform (argument)
	 * @return two matrix (first - with  amplitude, second - with phase)
	 */
	public ArrayList<double[][]> fft(double[][] amplMatrix, double[][] phaseMatrix){
		Compl[][] matrix = null;
		ArrayList<double[][]> result = new ArrayList<double[][]>();
		int height = amplMatrix[0].length;
		int width = amplMatrix.length;
		matrix = new Compl[width][height];
		for (int i = 0; i < width; i++) {
			for (int j = 0; j < height; j++) {
				matrix[i][j] = new Compl();
				matrix[i][j].setRe(amplMatrix[i][j]
						* Math.cos(phaseMatrix[i][j]));
				matrix[i][j].setIm(amplMatrix[i][j]
						* Math.sin(phaseMatrix[i][j]));
			}
		}
		fft(matrix); 
		double[][] resultAmpl = new double[width][height];
		double[][] resultPhase = new double[width][height];
		
		for (int i = 0; i < width; i++) {
			for (int j = 0; j < height; j++) {
				double module = Math.hypot(matrix[i][j].getRe(),
						matrix[i][j].getIm());// модуль
				double arg = Math.atan2(matrix[i][j].getRe(), matrix[i][j].getIm());
				resultAmpl[i][j] = module;
				resultPhase[i][j] = arg;
			}
		}
		result.add(0,resultAmpl);
		result.add(1,resultPhase);
		return result;
	}
	public ArrayList<BufferedImage> fft(AbstractDoe doe, DoeType type) {
		ArrayList<BufferedImage> imageList = new ArrayList<BufferedImage>();
		if (DoeType.AMPLITUDE == type) {
			amplMatrix = doe.getAbsoluteAmplitudeMatrix();
			phaseMatrix = doe.getSignPhaseMatrix();
		} else if (DoeType.PHASE == type) {
			phaseMatrix = doe.getPhaseTo2PIMatrix();
			amplMatrix = new double[phaseMatrix.length][phaseMatrix[0].length];
			for (int i = 0; i < phaseMatrix.length; i++) {
				for (int j = 0; j < phaseMatrix[0].length; j++) {
					amplMatrix[i][j] = 1;
				}
			}
		}
		else if (DoeType.PHASE_MODAN_AMPLITUDE == type) {
			amplMatrix = doe.getAbsoluteAmplitudeMatrix();
			phaseMatrix = doe.getPhaseTo2PIMatrix();
		}
		ArrayList<double[][]> afterFourer = fft(amplMatrix, phaseMatrix);
		resultAmpl = afterFourer.get(0);
		resultPhase = afterFourer.get(1);
		
		int height = resultAmpl[0].length;
		int width = resultAmpl.length;
		double maxAmpl = -1;
		double maxPhase = -1;
		BufferedImage imgAmplitude = new BufferedImage(width, height,
				BufferedImage.TYPE_INT_RGB);
		BufferedImage imgPhase = new BufferedImage(width, height,
				BufferedImage.TYPE_INT_RGB);
		for (int i = 0; i < width; i++) {
			for (int j = 0; j < height; j++) {
				maxAmpl = resultAmpl[i][j] > maxAmpl ? resultAmpl[i][j]
						: maxAmpl;
				maxPhase = resultPhase[i][j] > maxPhase ? resultPhase[i][j]
						: maxPhase;
			}
		}
		if (DoeType.AMPLITUDE == type) {
			for (int i = 0; i < width; i++) {
				for (int j = 0; j < height; j++) {
					resultAmpl[i][j] = resultAmpl[i][j] / maxAmpl * 255;// Integer.MAX_VALUE;
					imgAmplitude.setRGB(i, j, (int) (resultAmpl[i][j]));
					resultPhase[i][j] = resultPhase[i][j] > 0 ? 255 : 0f;
					imgPhase.setRGB(i, j, (int) (resultPhase[i][j]));
				}
			}
		}
		else if (DoeType.PHASE == type || DoeType.PHASE_MODAN_AMPLITUDE == type) {
			for (int i = 0; i < width; i++) {
				for (int j = 0; j < height; j++) {
					resultAmpl[i][j] = resultAmpl[i][j] / maxAmpl *  255;//Integer.MAX_VALUE;
					imgAmplitude.setRGB(i, j, (int) (resultAmpl[i][j]));
					resultPhase[i][j] = resultPhase[i][j] / maxPhase * 255;
					imgPhase.setRGB(i, j, (int) (resultPhase[i][j]));
				}
			}
		}
		imageList.add(0, imgAmplitude);
		imageList.add(1, imgPhase);
		return imageList;
	}
	public ArrayList<double[][]> fftMatrixes(AbstractDoe doe, DoeType type) {
		if (DoeType.AMPLITUDE == type) {
			amplMatrix = doe.getAbsoluteAmplitudeMatrix();
			phaseMatrix = doe.getSignPhaseMatrix();
		} else if (DoeType.PHASE == type) {
			phaseMatrix = doe.getPhaseTo2PIMatrix();
			amplMatrix = new double[phaseMatrix.length][phaseMatrix[0].length];
			for (int i = 0; i < phaseMatrix.length; i++) {
				for (int j = 0; j < phaseMatrix[0].length; j++) {
					amplMatrix[i][j] = 1;
				}
			}
		}
		else if (DoeType.PHASE_MODAN_AMPLITUDE == type) {
			amplMatrix = doe.getAbsoluteAmplitudeMatrix();
			phaseMatrix = doe.getPhaseTo2PIMatrix();
		}
		ArrayList<double[][]> afterFourer = fft(amplMatrix, phaseMatrix);
		resultAmpl = afterFourer.get(0);
		resultPhase = afterFourer.get(1);	
		return afterFourer;
	}

	/*
	 * public void fft(Compl[][] matrix, DoeType type) {
	 * 
	 * }
	 */

	private static Compl[][] fft(Compl[][] array)/*
												 * , DftNormalization norm,
												 * TransformType type)
												 */{
		int step1 = 0;
		int step2 = 0;
		int len1 = array.length;
		int len2 = array[0].length;
		if (len1 % 2 == 0 && len2 % 2 == 0) {
			while (len1 != 1) {
				if (len1 % 2 != 0) {
					throw new IllegalArgumentException(
							"Размер области не является степенью двойки");
				}
				len1 /= 2;
				step1++;

			}
			while (len2 != 1) {
				if (len2 % 2 != 0) {
					throw new IllegalArgumentException(
							"Размер области не является степенью двойки");
				}
				len2 /= 2;
				step2++;

			}
		}
		for (int i = 0; i < array.length; i++) {
			double[][] dataRI = new double[2][array[0].length];
			Compl[] temp = new Compl[array[0].length];
			for (int j = 0; j < array[0].length; j++) {
				temp[j] = array[i][j];
			}
			QuickFourierTransform.fftdhts(temp);
			for (int j = 0; j < array[0].length; j++) {
				array[i][j] = temp[j];
			}
		}

		for (int i = 0; i < array[0].length; i++) {
			double[][] dataRI = new double[2][array[0].length];
			Compl[] temp = new Compl[array.length];
			for (int j = 0; j < array.length; j++) {
				temp[j] = array[j][i];
			}
			QuickFourierTransform.fftdhts(temp);
			for (int j = 0; j < array.length; j++) {
				array[j][i] = temp[j];
			}
		}

		return inverseMatrix(array);
	}

	private static Compl[][] fft2D(Compl[][] array) {
		DoubleFFT_2D fft = new DoubleFFT_2D(array.length, array[0].length);
		double[][] a = new double[array.length][2 * array[0].length];
		for (int k1 = 0; k1 < array.length; k1++) {
			for (int k2 = 0; k2 < array[0].length; k2++) {
				a[k1][2 * k2] = array[k1][k2].getRe();
				a[k1][2 * k2 + 1] = array[k1][k2].getIm();
			}
		}
		fft.complexForward(a);
		for (int k1 = 0; k1 < array.length; k1++) {
			for (int k2 = 0; k2 < array[0].length; k2++) {
				array[k1][k2].setRe(a[k1][2 * k2]);
				array[k1][k2].setIm(a[k1][2 * k2 + 1]);
				// System.out.println(a[k1][2*k2] + "    " +a[k1][2*k2+1]);
			}
		}
		return array;
	}

	public static Compl[][] inverseMatrix(Compl[][] matrix) {
		Compl[][] result = new Compl[matrix.length][matrix[0].length];
		for (int k1 = 0; k1 < result.length / 2; k1++) {
			for (int k2 = 0; k2 < result[0].length / 2; k2++) {
				result[matrix.length / 2 + k1][matrix[0].length / 2 + k2] = matrix[k1][k2];
			}
			for (int k2 = result[0].length / 2; k2 < result[0].length; k2++) {
				result[k1 + matrix.length / 2][k2 - matrix[0].length / 2] = matrix[k1][k2];
			}
		}
		for (int k1 = result.length / 2; k1 < result.length; k1++) {
			for (int k2 = 0; k2 < result[0].length / 2; k2++) {
				result[k1 - matrix.length / 2][k2 + matrix[0].length / 2] = matrix[k1][k2];
			}
			for (int k2 = result[0].length / 2; k2 < result[0].length; k2++) {
				result[k1 - matrix.length / 2][k2 - matrix[0].length / 2] = matrix[k1][k2];
			}
		}
		return result;
	}

	private static Compl[] newFourer(Compl[] array, Integer nn,
			Boolean inverseFFT) {

		int ii, jj, n, mmax, m, j, istep, i, isign;
		double wtemp, wr, wpr, wpi, wi, theta, tempr, tempi;
		double[] a = new double[2 * array.length];
		for (int k2 = 0; k2 < array.length; k2++) {
			a[2 * k2] = array[k2].getRe();
			a[2 * k2 + 1] = array[k2].getIm();
		}

		if (inverseFFT) {
			isign = -1;
		}

		else {
			isign = 1;
		}

		n = 2 * nn;

		j = 1;

		ii = 1;

		while (ii <= nn) {

			i = 2 * ii - 1;

			if (j > i) {
				tempr = a[j - 1];

				tempi = a[j];

				a[j - 1] = a[i - 1];

				a[j] = a[i];

				a[i - 1] = tempr;

				a[i] = tempi;

			}

			m = (n - n % 2) / 2;

			while ((m >= 2) && (j > m)) {

				j = j - m;

				m = (m - m % 2) / 2;

			}

			j = j + m;

			ii++;
		}

		mmax = 2;

		while (n > mmax) {
			istep = 2 * mmax;

			theta = 2 * Math.PI / (isign * mmax);

			wpr = -2.0 * Math.sqrt(Math.sin(0.5 * theta));

			wpi = Math.sin(theta);

			wr = 1.0;

			wi = 0.0;

			ii = 1;

			while (ii <= (mmax - mmax % 2) / 2) {

				m = 2 * ii - 1;

				jj = 0;

				while (jj <= ((n - m) - (n - m) % istep)) {

					i = m + jj * istep;

					j = i + mmax;

					tempr = wr * a[j - 1] - wi * a[j];

					tempi = wr * a[j] + wi * a[j - 1];

					a[j - 1] = a[i - 1] - tempr;

					a[j] = a[i] - tempi;

					a[i - 1] = a[i - 1] + tempr;

					a[i] = a[i] + tempi;

					jj++;

				}

				wtemp = wr;

				wr = wr * wpr - wi * wpi + wr;

				wi = wi * wpr + wtemp * wpi + wi;

				ii++;

			}

			mmax = istep;

		}

		if (inverseFFT) {

			for (int I = 0; I < 2 * nn - 1; I++) {

				a[I - 1] = a[I - 1] / nn;

			}
		}
		for (int k2 = 0; k2 < array.length; k2++) {
			array[k2].setRe(a[2 * k2]);
			array[k2].setIm(a[2 * k2 + 1]);
		}
		return array;
	}

	private static Compl[] fftWiki(Compl[] input, boolean DIRECT) {
		// - n is the dimension of the problem
		// - nu is its logarithm in base e
		int n = input.length;

		// If n is a power of 2, then ld is an integer (_without_ decimals)
		double ld = Math.log(n) / Math.log(2.0);

		// Here I check if n is a power of 2. If exist decimals in ld, I quit
		// from the function returning null.
		if (((int) ld) - ld != 0) {
			System.out.println("The number of elements is not a power of 2.");
			return null;
		}

		// Declaration and initialization of the variables
		// ld should be an integer, actually, so I don't lose any information in
		// the cast
		int nu = (int) ld;
		int n2 = n / 2;
		int nu1 = nu - 1;
		double[] xReal = new double[n];
		double[] xImag = new double[n];
		double tReal, tImag, p, arg, c, s;

		// Here I check if I'm going to do the direct transform or the inverse
		// transform.
		double constant;
		if (DIRECT)
			constant = -2 * Math.PI;
		else
			constant = 2 * Math.PI;

		// I don't want to overwrite the input arrays, so here I copy them. This
		// choice adds \Theta(2n) to the complexity.
		for (int i = 0; i < n; i++) {
			xReal[i] = input[i].getRe();
			xImag[i] = input[i].getIm();
		}

		// First phase - calculation
		int k = 0;
		for (int l = 1; l <= nu; l++) {
			while (k < n) {
				for (int i = 1; i <= n2; i++) {
					p = bitreverseReference(k >> nu1, nu);
					// direct FFT or inverse FFT
					arg = constant * p / n;
					c = Math.cos(arg);
					s = Math.sin(arg);
					tReal = xReal[k + n2] * c + xImag[k + n2] * s;
					tImag = xImag[k + n2] * c - xReal[k + n2] * s;
					xReal[k + n2] = xReal[k] - tReal;
					xImag[k + n2] = xImag[k] - tImag;
					xReal[k] += tReal;
					xImag[k] += tImag;
					k++;
				}
				k += n2;
			}
			k = 0;
			nu1--;
			n2 /= 2;
		}

		// Second phase - recombination
		k = 0;
		int r;
		while (k < n) {
			r = bitreverseReference(k, nu);
			if (r > k) {
				tReal = xReal[k];
				tImag = xImag[k];
				xReal[k] = xReal[r];
				xImag[k] = xImag[r];
				xReal[r] = tReal;
				xImag[r] = tImag;
			}
			k++;
		}

		// Here I have to mix xReal and xImag to have an array (yes, it should
		// be possible to do this stuff in the earlier parts of the code, but
		// it's here to readibility).
		Compl[] newArray = new Compl[xReal.length];
		double radice = 1 / Math.sqrt(n);
		for (int i = 0; i < newArray.length; i++) {
			// int i2 = i / 2;
			// I used Stephen Wolfram's Mathematica as a reference so I'm going
			// to normalize the output while I'm copying the elements.
			newArray[i] = new Compl();
			newArray[i].setRe(xReal[i] * radice);
			newArray[i].setIm(xImag[i] * radice);
		}
		return newArray;
	}

	/**
	 * The reference bitreverse function.
	 */
	private static int bitreverseReference(int j, int nu) {
		int j2;
		int j1 = j;
		int k = 0;
		for (int i = 1; i <= nu; i++) {
			j2 = j1 / 2;
			k = 2 * k + j1 - 2 * j2;
			j1 = j2;
		}
		return k;
	}

	// *********************************************************************
	private static Compl[] fftIntrocs(Compl[] temp) {
		Complex[] x = new Complex[temp.length];
		for (int i = 0; i < temp.length; i++) {
			x[i] = new Complex(temp[i].getRe(), temp[i].getIm());
		}
		x = fftIntrocs(x);
		for (int i = 0; i < temp.length; i++) {
			temp[i] = new Compl(x[i].re(), x[i].im());
		}
		return temp;
	}

	private static Complex[] fftIntrocs(Complex[] x) {

		int N = x.length;

		// base case
		if (N == 1)
			return new Complex[] { x[0] };

		// radix 2 Cooley-Tukey FFT
		if (N % 2 != 0) {
			throw new RuntimeException("N is not a power of 2");
		}

		// fft of even terms
		Complex[] even = new Complex[N / 2];
		for (int k = 0; k < N / 2; k++) {
			even[k] = x[2 * k];
		}
		Complex[] q = fftIntrocs(even);

		// fft of odd terms
		Complex[] odd = even; // reuse the array
		for (int k = 0; k < N / 2; k++) {
			odd[k] = x[2 * k + 1];
		}
		Complex[] r = fftIntrocs(odd);

		// combine
		Complex[] y = new Complex[N];
		for (int k = 0; k < N / 2; k++) {
			double kth = -2 * k * Math.PI / N;
			Complex wk = new Complex(Math.cos(kth), Math.sin(kth));
			y[k] = q[k].plus(wk.times(r[k]));
			y[k + N / 2] = q[k].minus(wk.times(r[k]));
		}
		return y;
	}

	// ******************************************

	private static Compl[] computeDft(Compl[] temp) {
		double[] inreal = new double[temp.length];
		double[] inimag = new double[temp.length];
		double[] outreal = new double[temp.length];
		double[] outimag = new double[temp.length];
		for (int i = 0; i < temp.length; i++) {
			inreal[i] = temp[i].getRe();
			inimag[i] = temp[i].getIm();
		}
		computeDft(inreal, inimag, outreal, outimag);
		for (int i = 0; i < temp.length; i++) {
			temp[i] = new Compl(outreal[i], outimag[i]);
		}
		return temp;
	}

	// http://nayuki.eigenstate.org/page/how-to-implement-the-discrete-fourier-transform
	private static void computeDft(double[] inreal, double[] inimag,
			double[] outreal, double[] outimag) {
		int n = inreal.length;
		// System.out.println("Start");
		for (int k = 0; k < n; k++) { // For each output element
			double sumreal = 0;
			double sumimag = 0;
			for (int t = 0; t < n; t++) { // For each input element
				sumreal += inreal[t] * Math.cos(2 * Math.PI * t * k / n)
						+ inimag[t] * Math.sin(2 * Math.PI * t * k / n);
				sumimag += -inreal[t] * Math.sin(2 * Math.PI * t * k / n)
						+ inimag[t] * Math.cos(2 * Math.PI * t * k / n);
			}
			outreal[k] = sumreal;
			outimag[k] = sumimag;
			// System.out.println(k);
		}
		// System.out.println("End");
	}
	public double[][] getAmplMatrix() {
		return amplMatrix;
	}
	public double[][] getPhaseMatrix() {
		return phaseMatrix;
	}
	public double[][] getResultAmpl() {
		return resultAmpl;
	}
	public double[][] getResultInt() {
		double[][] res =  resultAmpl;
		for(int i =0; i<res.length;i++)
		{
			for(int j=0;j<res[0].length;j++)
			{
				res[i][j]=res[i][j]*res[i][j];
			}
		}
		return res;
	}
	public double[][] getResultPhase() {
		return resultPhase;
	}

}